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Structural and transport properties of interacting localized flux lines in the Bose glass phase of 
irradiated superconductors are studied by means of Monte Carlo simulations near the matching 
field B§, where the densities of vortices and columnar defects are equal. For a completely random 
columnar pin distribution in the icy-plane transverse to the magnetic field, our results show that 
the repulsive vortex interactions destroy the Mott insulator phase which was predicted to occur at 
B — 5$. On the other hand, for ratios of the penetration depth to average defect distance X/d < 1, 
characteristic remnants of the Mott insulator singularities remain visible in experimentally accessible 
quantities as the magnetization, the bulk modulus, and the magnetization relaxation, when B is 
varied near B§. For spatially more regular disorder, e.g., a nearly triangular defect distribution, we 
find that the Mott insulator phase can survive up to considerably large interaction range X/d, and 
may thus be observable in experiments. 

PACS numbers: 74.60. Ge, 05.60,+w 



I. INTRODUCTION 

The remarkably rich phase diagram of magnetic flux 
lines in high-T c superconductors, especially when sub- 
ject to point and/or extended disorder, has attracted 
considerable experimental and theoretical interest .u For 
the purpose of applying type-II superconductors in ex- 
ternal magnetic fields, an effective vortex pinning mech- 
anism is essential, in order to minimize dissipative losses 
caused by the Lorentz-force induced movement of flux 
lines across the sample. This issue becomes even more 
important for the high-T c superconductors, because the 
strongly enhanced thermal fluctuations in these materi- 
als render the Abrikosov flux lattice unstable. Thus, well 
below the upper critical field H C2 (T), there occurs a first- 
order melting transition leading to a normal-conducting 
flux liquid phase. □ It turns out to be rather difficult to 
pin such a vortex liquid consisting of strongly fluctuat- 
ing flexible lines by just intrinsic point disorder (e.g., 
oxygen vacancies), although asymptotically a truly su- 
perconducting low-temperature vortex glass phase was 
predicted. u 

Therefore, in addition to point defects, the influence of 
extended or correlated disorder, promising stronger pin- 
ning effects, on vortex transport properties has been con- 
sidered. Experimentally, linear damage tracks have been 
produced in oxide superconductors by irradiation with 
high-energy heavy ions. These columnar defects serve 
as strong pinning centers for the flux lines, by signifi- 
cantly increasing the critical current and shifting the ir- 
reversibility line upwards (see Ref. [ [jj , Sec. IX B for refer- 
ences). Theoretically, the statistical mechanics problem 
of directed lines interacting with columnar pins can be 
mapped onto the quantum mechanics of two-dimensional 
bosons subject to static point disorder via the identifica- 
tion of the vortex trajectories r^(z) with particle world 



lines in imaginary timeau □. In this mathematical anal- 
ogy, thermal fluctuations ~ fceT map onto quantum fluc- 
tuations ~ h, the vortex line tension e\ represents the 
boson mass, and the sample thickness L along the z di- 
rection parallel to the magnetic field B corresponds to 
the inverse temperature /3h of the quantum system (such 
that the thermodynamic limit L — > oo implies zero tem- 
perature for the fictitious bosons) . The boson picture can 
then be employed to immediately construct the expected 
phase diagram for flux lines interacting with columnar 
defects, as depicted in Fig. [TP'S and its properties can 
be analyzed by utilizing the (sgaling) theory of bosons 
subject to static point disorder .0 

At high temperatures (and ia the thermodynamic limit 
L — > oo, i.e., for thick samplesu), one thus finds an entan- 
gled liquid of unbound flux lines (corresponding to a bo- 
son superfluid), separated by a sharp second-order tran- 
sition from a low-temperature phase of localized vortices. 
This Bose glass phase is characterized by an infinite tilt 
modulus C44, and turns out to be stable over a certain 
range of tipping angles of B away from the z direction 
(transverse Meissner effect). □ At low temperatures in the 
Bose glass, the dominant vortex transport mechanism 
for low external currents is expected to be the analog 
of variable-range hopping in disordered semiconductors!!! 
namely the formation of double kinks extending from one 
columar defect to another energetically favorable, but 
not necessarily adjacent pinning site.El Therefore the ac- 
tual distribution of pinning energies constitutes an im- 
portant feature of the system that determines its low- 
energy current-voltage (I-V) characteristics. As in the 
case of doped semiconductors p~t3 long-range repulsive 
interactions drastically modify the distribution of pin- 
ning energies, which is the analog of an (interacting) 
density of states, and may lead to the emergence of a 
correlation-induced soft "Coulomb gap" separating the 
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FIG. 1. Phenomenological phase diagram of a type-II su- 
perconductor with strong correlated disorder. A sharp phase 
transition line T C (B) separates the flux liquid from the Bose 
glass phase. Within this phase, one may identify a crossover 
line discriminating between the weakly and strongly pinned 
Bose glass. The Mott insulator is located deep within the 
Bose glass at low temperatures, appearing as a line phase 
precisely at the matching field B = _B$. 

filled and empty states .El This effect in turn reduces vor- 
tex transport in the variable-range hopping regime con- 
siderably. Recent experiments have indeed identified the 
characteristic signatures of variable-xangp hopping pro- 
cesses in irradiated superconductors.EiJO 

At least for weakly interacting flux lines, the theory 
also suggests a low-temperature Mott insulator phase, 
occurring when B = B$, i.e., the vortex density exactly 
matches that of the columnar damage tracks!] In this 
situation, all the defects will be occupied by one vortex 
each, and the introduction of another flux line into the 
system will require the typical pinning energy Uq (pro- 
vided the repulsive vortex interactions are negligible — 
weak repulsive forces will slightly increase this energy 
cost). Consequently, upon increasing the external field 
H, one expects a "lock- in" effect at B — B$, where the 
internal field remains constant until this energy barrier is 
overcome, as shown in Fig. |(seealsoRefs. [0^o|). The 
Mott insulator phase is therefore characterized by a hard 
gap in the distribution of pinning energies B Correspond- 
ingly, the (reversible) magnetization M = B — H should 
display a sharp downward jump exactly at B = this 
can also be seen from the relation M = —dG/dB, where 
G is the total Gibbs free energy of the system, here domi- 
nated by the internal energy E, which increases abruptly 
at £?$ as a new vortex enters the sample. The Bose 
glass and the Mott insulator represent distinct thermo- 
dynamic phases, for the latter should be characterized 
by an infinite tilt modulus C44, and a diverging compres- 
sional modulus en ~ dH/dB. Characteristic minima in 
the reversible magnetization near the matching field were 
indeed observed experimentally at low temperatures 
single-crystal Tl coBapounds,EZl and both BSCCO tapesl 
and single-crystals.E3 Furthermore, pronounced minima 
in the magnetization relaxation at B w 1.4B$ were re- 
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FIG. 2. Left: B vs H showing the typical Mott insulator 
"lock-in" at the matching field B§. B does not change over a 
certain range of the applied field [H a , Hb] , until the magnetic 
field energy (at Hb) becomes large enough to push another 
vortex into the system. Right: At the Mott insulator transi- 
tion, the magnetization M = B — H should display a sharp 
jump of size oc (Hb — H a ) as a result of this "lock-in" effect. 



ported for YBCO single crystals ,E3-and for B m I.4£?$ 
and B w 3.3_B$ in a Tl compound)^ and interpreted as 
further evidence for vortex "lock-in" behavior character- 
istic of the Mott insulator phase. 

One of the main results of our numerical study will be, 
however, that the original predictions of a distinct ther- 
modynamic Mott insulator phase are, strictly speaking, 
applicable only in two limiting cases: (i) For short-range 
interactions, i.e., when the London penetration depth A 
is considerably smaller than the average defect spacing d, 
and the sole effect of the intervortex repulsion is to pre- 
vent double occupancy of columnar pins; but it cannot 
prevent all the defect sites, even for a spatially random 
distribution of pinning centers, to be filled with vortices, 
(ii) For a regular array of columnar defects, e.g., a trian- 
gular or square lattice, for which because of the spatial 
symmetry the interactions merely lead to a constant ad- 
ditive renormalization of the pinning energies, equal for 
all defect lattice sites. In all intermediate cases, i.e., a 
disordered spatial defect distribution and a sizeable in- 
teraction range X/d, one has to expect that it will be im- 
possible to fill the columnar pins completely. For, the oc- 
cupation of very close pins will be energetically costly as 
compared to moving a vortex into a high-symmetry "in- 
terstitial" position, at maximum distance from neighbor- 
ing defect sites. Thus, a departure from both situations 
(i) and (ii) by either increasing the interaction strength 
and range, or introducing disorder into the defect con- 
figuration, will lead to a softening of the ideal Mott in- 
sulator hard gap in the distribution of pinning energies, 
and smear the corresponding singularities in the free en- 
ergy G, magnetization M, and bulk modulus c\\. In the 
thermodynamic sense, the Mott insulator does therefore 
not exist as a distinct phase in typical experimental sit- 
uations. On the other hand, if X/d is not too large, we 
shall find that very characteristic remnants of the Mott 
insulator singularities appear in measurable quantities, 
and we shall indeed be able to explain why and in which 
direction these maxima or minima are shifted away from 
B = _B$, as experimentally observed.t£l Also, recently 
regular arrays of columnar defects have been produced 
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in high-T c materialsj=9 and in these more artificial sam- 
ples the Mott insulator phase might play a much more 
prominent role. 

Another important effect of the vortex interactions 
leads to the distinction of a strongly and a weakly pinned 
regime, respectively, within the Bose glass phase. QiiS'EZl 
For short-range interactions X/d <C 1, or for larger X/d, 
but low filling / = B/Bq, < 0.5, the vortices will be 
strongly bound to the columnar defects, in the latter 
situation because of the collective "Coulomb gap" cor- 
relation effects which render the transfer of a flux-Line 
to another equally favorable site extremely unlikely.EJ In 
the other extreme case, namely at large fillings / > I and 
for very strong interactions, the pinning potentials actu- 
ally constitute a rather weak perturbation of the strongly 
correlated vortex system (with the flux lines arranging 
themselves as closely as possible to a triangular lattice), 
and the entire system can be shifted easily through the 
application of an external current. However, there should 
also be an intermediate weakly pinned Bose glass regime 
around the matching field, where some of the flux lines 
are bound to the strong pinning sites, and the others are 
held in place through tJaeic, mutual repulsive interactions 
with adjacent vortices£3'E3 In a sample with a spatially 
random distribution of columnar defects at low temper- 
atures, one would therefore expect a crossover from the 
Mott insulator to a weakly pinned Bose glass, as X/d is 
increased. Estimates for the location of this crossover 
line (denoted B* (T) there) are given in Ref . [ [| . 

Thus, if one aims to design a high-T c superconduct- 
ing sample in such a way that its pinning properties are 
optimized, clearly one needs to understand the influence 
of the repulsive interactions and the effect of different 
spatial pinning site distributions in a quantitative man- 
ner. Obviously, there are two promising regimes for such 
a specific material "tailoring" , namely either the under- 
filled (/ < 0.5) strongly pinned Bose glass where the en- 
suing soft correlation gap in the distribution of pinning 
energies ensures that at least in the variable-range hop- 
ping regime vortex transport is effectively reduced, or the 
incompressible Mott insulator phase with its hard energy 
barrier preventing flux line motion. It is the purpose of 
this paper, therefore, to provide a quantitative numer- 
ical study of the low-temperature structural and trans- 
port properties of the entire Bose glass phase, including 
and specifically addressing the region near the matching 
field B = B&, as well as the features of the intermediate 
weakly pinned regime. Our investigations are based on 
a variant of the Coulomb glass model, which we modify 
in such a way that high-symmetry interstitial sites are 
accessible to the vortices in addition to the defect posi- 
tions themselves, and we may thus utilize the successful 
algorithms and. considerable experience obtained in ear- 
lier studies.ErEj This model will be described in detail in 
Sec. H while in Sec. |ffl we shall introduce and define the 
quantities of interest to characterize both relevant static 
and dynamic properties of the Bose glass and Mott in- 
sulator phases in irradiated superconductors. Sec. IV is 



contains the main results of our numerical study. We 
shall explore both a random distribution of pinning sites 
(uniform disorder) and spatially correlated disorder dis- 
tributions, as obtained by deformations of a perfect tri- 
angular defect lattice. Finally, we shall summarize and 
discuss our results. 



II. THE MODEL 

The theoretical modeling of the Bose glass phase is 
based on the following free energy for Ny flux lines, 
described by their two-dimensional trajectories r^(z) as 
they traverse the sample of thickness L, interacting with 
each other and with Ng, columnar defects parallel to the 
magnetic field B || z,P 



r L Nv \ ~ / 7 / \ \ 2 , N v 



£Vb[ri(z)-R fc ] 



(1) 



The first term in Eq. ([!]) describes the elastic flux line 
tension and originates in an expansion of the line energy 
of nearly straight vortices with respect to small tipping 
angles. □ The tilt modulus is given by e\ ~ eo ln(ao/£)j 
where the energy scale is set by e = (^o/^ttA) 2 , a = 
(4/3)V4(0 Q /s) i s the average lattice spacing between the 
flux lines, £ is the superconducting correlation length, 
and A denotes the London penetration depth. 

The second term represents the interaction of all vortex 
pairs (which we approximate to be local in z), where ry = 
|rj — Yj\ and V(r) = 2eoK (r/X) is the screened repul- 
sive vortex potential, with the modified Bessel function 
Kq(x) ~ — ln(x) as x — > 0, and K$(x) ~ x" 1 / 2 exp(— x) 
for x — > oo. 

Finally, the No columnar pins || B are modeled 
by a sum of z-independent attractive potential wells 
Vjj(r — Rfc) = — Uk 9(|r — Rfc|), with average spacing 
d, and centered at randomly distributed positions {RaJ. 
Here, O denotes the Heaviside step function. The damage 
track diameters, produced by heavy ion irradiation, are 
typically 2co ~ 100 A wide, with a variation induced by 
the root-mean-square energy dispersion of the ion beam 
of Sck/co « 15%, where c k denotes the individual defect 
radii. The pinning energies Uk are related to the column 
diameters via. the interpolation formula given by Nelson 
and VinokurD 



C/ fc «|ln[l + (c fc /V2e) 2 ] 



(2) 



where £ w 10 A is the coherence length. The distribu- 
tion of column diameters thus induces a variation of the 
pinning energies with the width 



then devoted to the analysis of our simulation data and 
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Above a certain temperature To, the effective radius 
of the normal-conducting region at the defects will be 
given by the coherence length, which in mean-field the- 
ory grows with temperature according to 



f(T) = Co(l-T/T £ 



x-l/2 



Hence, y/2£(T) should be inserted into Eq. (J 
Cfc, whenever y/2£(T) > Co, or T > To, with 

T = T c (l - 2$/c§) . 



(4) 



instead of 



(5) 



Consequently, above To the fluctuations of the pinning 
energies will be effectively smoothed out. With £ ~ 10 
A and c w 50 A this will happen at T sa 0.92T C . As we 
shall describe later, we will not be particularly interested 
in this regime, since the Mott insulator is predicted to 
exist deep in the Bose glass phase at low temperatures 
when thermal fluctuations can be neglected. 

The classical statistical mechanics properties of Eq. (|l|) 
can be derived from an analysis of the grand-canonical 
partition functionQ (we set the Boltzmann constant fee 
1 throughout this paper) 



oo 

Z = E ^exp( M JV K /T)x 

JVv=0 V ' 
, N v 

J ndr«( Z )e X p{-^[r,(z)]/T} 



where fi denotes the chemical potential 



<t>0 „ i ^0 



(H-H C1 ) , 



(6) 



(7) 



which changes sign at the lower critical field H Cl = 
<f>o In k/(47tA 2 ). The analytic approach to this parti- 
tion function has made use of a mapping onto a two- 
dimensional quantum-mechanical problem with a static 
(i.e., z-independent) disorder potential, where the physi 
cal direction z has been interpreted as imaginary time 
In the thermodynamic limit L — > oo (which tunes the ef- 
fective temperature of the Bose system to zero), the prop- 
erties of the model are determined by the lowest-energy 
eigenstate of the corresponding transfer matrix or equiv- 
alent quantum-mechanical problem. The corresponding 
ground state wave function is symmetric with respect to 
exchange of flux lines, and thus of bosonic character. The 
ensuing localization transition found with this method, 
where the bosons become bound to the static disorder 
potential (i.e., the flux lines are pinned to the columnar 
defects), leads to a disorder-dominated low-temperature 
phase, termed Bose glass. 

Here, we will be interested in the low-temperature 
(ground state, T = 0) properties of Eq. ([!]), which is 
a fair approximation to the Bose glass phase for temper- 
atures T < T\ . Ti denotes the temperature above which 
entropic corrections associated with thermal fluctuations 
become relevant for the pinned flux lines. Physically, this 



means that vortices start to fluctuate around their de- 
fect positions and can no longer be thought of as straj; 
lines. Theoretical estimates yield T\ ~ 0.6 . . . 0.8T C . 
although in some expe riments T\ has been found to be 
considerably lower .PtM At an even higher depinning tem- 
perature T\ < Tdp < T c , vortices start to leave the 
columnar defects, but may still be effectively pinned 
by their mutual interactions with the neighboring occu- 
pied defects. We restrict our study to the regime be- 
low Ti, where the vortices are essentially straight lines, 
which simplifies our problem considerably. When com- 
paring our simulation results later with experimental 
data, T < 0.6 . . . 0.8T C will a posteriori prove to be a 
good approximation. 

Since thermal wandering and bending of vortex lines 
below T\ is suppressed, we can neglect the bending en- 
ergy term in Eq. ([!]). This leaves us with a time- 
independent, two-dimensional problem of Ny interact- 
ing "particles" and No defects. We also aim to address 
a possible Mott insulator phase within the Bose glass, 
near the matching field _B$. Consequently, we want to 
study a filling fraction 



Ny 



B 



1 



(8) 



which is why we have to account for the possibility (es- 
pecially for / > 1) that not all vortices can be accommo- 
dated on defect sites. Hence, we represent our problem 
on an underlying triangular grid with N "lattice" sites 
and N — Nu interdefect positions ( "interstitials" ) . This 
is the main new aspect of our simulations in compari- 
son to those reported in Ref. [ [l2|], where only defect 
sites were available. Thus, that study was limited to low 
filling fractions / <C 1. The effective, two-dimensional 
Hamiltonian used in our simulations, therefore reads 



1 N N D 



(9) 



fc=l 



and its grand-canonical counterpart is 



H 



eff 



mE'' 



(10) 



Here, n.j = {0, 1} represents the site occupations num- 
ber; note that 52»=i n i ~ ^v- Taking £ 10 A (at 
T = 0) and Co ~ 50 A, we obtain for the bare pinning 
energies i/. = —{Uk) + Wk, with (Uk) = 0.65 and width 
w = 0.1, both in units of 2eo- For simplicity, we assume 
a flat distribution of pinning energies around its mean, 
P{wk) = ®(w — \wk\)/(2w). Our choice of an underly- 
ing triangular grid is motivated by the fact that without 
disorder the energetically most favorable state for the 
vortices is the (triangular) Abrikosov lattice. Thus in 
effect we take into account the energetically lowest high- 
symmetry "interstitial" positions between the columnar 
defects. 
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Finding the ground state configuration of this sys- 
tem by analytical means constitutes a very difficult task, 
which has to date not been successfully mastered. Conse- 
quently, we have to resort to numerical simulation tech- 
niques with a suitable ground state finding algorithm. 
Fortunately, the Hamiltonian (g) is precisely of the form 
studied in the context of charge carriers localized at ran- 
dom impurities in doped semiconductors (i.e., the so- 
called Coulomb glass problem with V(r) ~ 1/r instead 
of our logarithmic repulsion), and we can utilize the 
broad numerical-experience gathered in those previous 
investigations B~V3 We shall use Eq. (g) in all of our sim- 
ulations, and will now introduce the simulation technique 
and the quantities we intend to measure. 



III. SIMULATION TECHNIQUE AND 
QUANTITIES OF INTEREST 

In this section we shall first describe in detail the 
Monte Carlo algorithm we have employed in all simu- 
lations presented in this paper. We then proceed with an 
introduction to the relevant quantities we have measured, 
and which can be qualitatively and semi-quantitatively 
compared with experimental data that aimed to identify 
a Mott insulator transition in YBCO, BSCCO, and Tl- 
compounds. 



A. Energy Minimization Procedure 

Our starting Hamiltonian is given in Eq. (^). We will 
be interested in the ground state properties of this Hamil- 
tonian as a function of the filling fraction /, which we 
shall vary near the predicted Mott insulator value / = 1 , 
and as a function of the interaction length scale A or, 
more often, as function of the dimensionless ratio X/d 
where d is the average defect distance. In order to find 
the ground state configuration of this Hamiltonian, we 
follow closely an energy minimization scheme described 
in detail by Shklovskii and Efros in Ref. [ 0], Chap. 14, 
with some-slight modifications according to Davies, Lee, 
and RiceEj Taking a triangular lattice with lattice con- 
stant a = 1 and N sites, and imposing periodic boundary 
conditions, we initially distribute Njj defects randomly 
on this underlying triangular grid. Each defect site is 
assigned a randomly chosen energy tk, from a flat dis- 
tribution with mean — (C/fc)/(2eo) = —0.65 and width 
u>/(2eo) = 0.1. Most of our simulations used N = 1600 
and No = 100, such that the average defect distance was 
d — 4, in units of the lattice constant a = 1. 

We then distribute Ny = fNo vortices randomly on 
the grid and follow the algorithm by Shklovskii and Efros 
in order to minimize the energy. This is done as follows: 
First, for our interacting system, we define single-particle 
site energies as 



dHcff 



N 



For a filled site, denotes the energy to remove a par- 
ticle from site i to infinity, while, for an empty site, 
is the energy needed to introduce an additional parti- 
cle from infinity on site i. In thermal equilibrium, the 
chemical potential fi energetically separates the occupied 
from the unoccupied states, i.e., q < fiii with rii = 1 
and ti > /iVi with = 0. In the random initial con- 
figuration, which usually corresponds to a strongly non- 
equilibrium situation, many occupied sites will have a 
higher energy then many unoccupied sites. To approxi- 
mate the ground state as closely as possible, we proceed 
in two steps. After the calculation of all we try to 
relax the initial state by successive particle-hole transi- 
tions. To do so, we determine the occupied site p with 
highest energy, e p = max.r„. =1 }. and the empty site q 
with lowest energy, e q = min{„. =0 } e i- If < £p> the 
particle on site p is moved to site q, the site energies are 
recalculated and the procedure is repeated until eventu- 
ally £p < e q \/p £ {i\rii = 1} , Vg € {i\rii — 0}. Thus we 
have now separated the occupied from the unoccupied 
states energetically and can give a first estimate of the 
chemical potential by 



(a + 



(12) 



The ensuing state is, however, still a very poor approx- 
imation of the ground state, as it is generally unstable 
against single-particle hops. To test this, we next com- 
pute the hopping energies 



V{n ) 



(13) 



This expression can be understood as follows. First, re- 
move a particle from site i to infinity, which yields an 
energy gain of ej. Then, take the particle back to site j, 
which costs the amount of tj — V{rij). The additional 
contribution of V(rij) stems from the fact that after re- 
moving the particle from site i there are only Ny — 1 
particles left, but ej is defined in terms of a system with 
Ny particles. Thus, the fictitious interaction with site i 
has to be accounted for explicitly. Note that a necessary 
condition for the ground state is A$_>j > for all pairs 
with rii = 1, rij = 0. 
Thus, in a second step we compute all Ai—>j, and then 
search for the minimum of all the negative hopping en- 



ergies: A p _ 



1{A^<0} 



Then, we perform 



fijV(rij) + U 



(11) 



the particle transfer from site p to site q, which reduces 
the total energy E of the system further. As this hop 
typically leads to a non-equilibrium state again, we next 
recalculate the site energies and repeat the particle-hole 
transfer procedure of the preceding paragraph. Then, 
another hopping transfer is attempted until eventually 
A,_>j > OVi 7^ j. The ensuing configuration is accepted 
as a fair approximation of the ground state for a partic- 
ular defect distribution, filling fraction / and ratio X/d. 
A typical initial and final particle distribution obtained 
with this procedure is shown in Fig. [| From the re- 
sulting configuration the site energies e^, the chemical 
potential /i, and the total energy E of the system can be 
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B. Distribution of (Interacting) Pinning Energies 



FIG. 3. Left: Typical initial configuration of randomly dis- 
tributed defects (o) and vortices (x) on a triangular lattice 
(dots) with N = 1600 and N v = N D = 100, i.e., / = 1; 
here we set X/d = 1. Right: Final configuration after the 
relaxation procedure. Most of the vortices are now located 
on energetically favorable defect sites, except for five parti- 
cles, since the respective defect sites are too close together, 
such that one of two/three close vortices is pushed away into 
low-energy "interstitial" voids due to the repulsive interaction 
V{n 3 ). 



computed. This procedure is repeated for a number of k 
different initial realizations of the disorder, such that we 
can obtain an ensemble average of the chemical potential 
\p], [E] and other quantities of interest, which we shall 
describe in the following subsections. We typically take 
k = 50 which has proven to be sufficient in our simula- 
tions. 

Of course, with this approximation scheme we still may 
not have reached the true ground state configuration, but 
rather a metastable "pseudo" -ground state. In princi- 
ple, one would have to test that pseudo-ground state 
against m = 2,3,...Ny particle hops, in order to ob- 
tain the true minimum-energy configuration. Previous 
investigations have shown, though, that terminating at 
m = 1 already .gives reliable estimates for the energy-level 
distribution.Erty To test the algorithm, we have equi- 
librated different initial configurations of vortices with 
an identical defect distribution. For example, we chose 
initially a) a completely random configuration, b) a tri- 
angular array, and c) we first occupied all defects with 
vortices, and then distributed the remaining vortices (if 
/ > 1) randomly, for given / and X/d. All different con- 
figurations yielded the same energies and fj, within the 
numerical uncertainties, although the spatial configura- 
tions typically looked slightly different. Additional tests 



of the algorithm will be presented in Sec. IV 



Having introduced our simulation algorithm, we will 
employ this —technique for 0.3 < / < 3 and X/d = 
1/4, 1/2, 1, 5c3in order to investigate whether there exists 
a Mott insulator within the Bosc glass phase, and see if 
it is stable against increasing the vortex- vortex repulsion 
as X/d is enhanced, while keeping the pinning strengths 
Uk fixed. In order to identify the Mott insulator and to 
compare with experimental data, we need to introduce 
a number of quantities that we can compute from our 
ground state configurations. This will be the topic of the 
next four subsections. 



From the interacting single-site energies defined in 
Eq. ([ll]) we can construct a histogram of all energies of 
the final ground state configuration. The quantity g(e)de 
then denotes the number of states per unit area with 
energies in the interval [e, e + de] with the normalizing 
constraint 



g{e)de = 1/a 2 



1 



(14) 



The quantity g(e) thus defines the interacting boson 
single-particle density of states or, in the flux line con- 
text, the distribution of interacting pinning energies. 
This quantity is particularly interesting in order to iden- 
tify a Mott insulator. For a true Mott insulator phase, 
5(e) should display a hard gap A near the chemical po- 
tential (j,, i.e., the (lower) energies of the occupied states 
should be separated from the (higher) energies of the un- 
occupied "band" by a region of finite width which con- 
tains no states at all. In the strong screening limit A — > 0, 
when there is no interaction between the vortices, the gap 
should be exactly A = (Uk) — w (see the top left of Fig. 
below) . 

Moreover, by further exploiting the analogy to two- 
dimensional bosons localized at randomly distributed de- 
fect sites, one has to consider the possibility that these 
spatial correlations produce a soft Coulomb gap in the 
distribution of pinning energies near the chemical po- 
tential fi. Such a Coulomb gap has indeed been found 
in the Coulomb glass problem of disordered semiconduc- 
tors, where the interactipiLis also long-range (~ 1/r up to 
the screening length A) J£tEJ and in the Bose glass phase 
for small filliags / and large X/d, i.e., almost logarithmic 



interactions Bt5 In both situations the density of states 
vanishes precisely at the chemical potential, and in its 
vicinity increases according to a power law, 



9(e) 



e-ii 



(15) 



Here, we have defined the effective gap exponent s c g, 
which depends strongly on A and on the filling fraction 
/, as we shall see. In fact, for small fillings / < 1 and 
short-range interactions A — * 0, the exponent is trivially, 
s c ff = 0, i.e., the density of states near fj, is constant J§tHl 
For A — > 00, the theoretical prediction within the frame- 
work of a mean-field type of analysis would be s e g = 00. 
With a random distribution of pinning sites, however, 
one has to expect that local energy fluctuations will play 
an important role, and thus < s c ff < 00 , with values 
that may-.depend on the filling /, and the actual energy 
scaleBH 

For finite X/d, and sufficiently low energies in the ut- 
most vicinity of //, one must find a crossover to s e g = 
in an infinite system. In a finite system, however, or at 
some distance from the chemical potential, one will find 
a nonzero effective exponent s e g, which can be obtained 
simply by fitting a power law in the vicinity of n (see also 
Rcf. [ Q). The determination of s e ff depending on X/d 



G 



and / is one of the tasks of our work since it affects cru- 
cially the experimentally measurable I-V characteristics 
of superconductors in the Bose glass phase, as we shall 
see below. Primarily, though, we want to see whether a 
hard gap A in the density of states survives the intro- 
duction of long-range interactions. 



Thus, the maximum of c\\ occurs at the location of the 
steepest negative slope in M(B), which usually is at a 
lower B value than the minimum in M. Only in the true 
Mott insulator phase will the singularities in en and M 
coincide precisely at B = 



C. Bulk Modulus 

From our measurement of the chemical potential p as 
a function of filling fraction / and constant interaction 
scale X/d, we can determine the bulk modulus en via 



en 



dp dp 
df- ND M^ 



(16) 



where k denotes the compressibility of the system. Phys- 
ically, the bulk modulus characterizes the stiffness of the 
system with respect to the insertion of an additional par- 
ticle into it. A divergence of en at the matching field 
B = Bg, (/ = 1) would indicate the presence of the Mott 
insulator phase, since the chemical potential p should 
show a sharp jump when one more vortex is put into 
the system. In the case of A — > 0, the jump becomes 
of the size of the hard gap in g(e), A = ([/&) — u>, i.e., 
introducing another vortex into the system costs a finite 
energy A. In this respect, the system becomes similar to 
a Meissner phase which has zero compressibility as well. 



D. Reversible Magnetization 

Since near the Mott insulator the Bose glass should 
become more and more diamagnetic due to the diverging 
bulk modulus, similar to the situation in the Meissner 
phase, a number of experiments have tried to identify 
the Mott insulator by measuring the magnetiz at i pB j-as a 
function of the (internal) magnetic field M(B)JlZE3E3 In 
our simulations we can also obtain the reversible mag- 
netization as a function of / from the total free energy 
of the system G, which becomes identical to the inter- 
nal energy E at T — 0, by noting the thermodynamic 
relation 



~ ~dB ~ ~AV W 



(17) 



Whereas in unirradiated samples the magneti: 
crease with / should be proportional to In 
the signature of a true Mott insulator would be a sharp 
negative dip in the magnetization at / = 1 due to the 
diamagnetic effect (diverging bulk modulus), see Fig. || 
above. Consequently, the bulk modulus and the magne- 
tization are in fact closely related to each other. Noting 
that M = B - H, / ~ B, and p ~ H - 



H Cl , we can write 



Cll 



dp 
df 



OH 
~dB 



dM 
~dB 



d 2 G 
dB 2 



(18) 



E. I-V Characteristics and Magnetization Relaxation 

The I-V characteristics measured in the Bose glass 
phase below T c have the same form as predicted for the 
vortex glass, namely 

E = po J cxp[-[/( J)/T] = p Ja q y[-{Uc/T){J 1 /JY^] , 

(19) 

where po is the normal state resistivity, U c is an en- 
ergy barrier scale, J\ a current density scale and most 
importantly, p e g denotes the effective glass or transport 
exponentn Experimentally, one determines the glass ex- 
ponent in different filling regimes and for different tem- 
peratures (effective interaction length scales A). We will 
see in this section that at low currents p c s is closely con- 
nected to s e ff (in the variable-range hopping regime, see 
below) via the (approximate) equation 



Pett 



Seff + 1 

D + s eS + l ' 



(20) 



where D is the dimension perpendicular to the applied 
field (here D = 2). In order to understand how we can 
extract the glass exponent from our data of the density of 
states, we first have to distinguish between two current 
density regimes. The discussion presented here essen- 
tially follows Refs. 



1. Hopping via vortex half-loop formation 

Let us first consider fairly strong currents J\ < J < J c , 
where J c is the critical depinning current above which 
superconductivity is destroyed. We shall see that in this 
regime, the motion of a single wortex is unaffected by the 
other flux lines in the sample A current perpendicular 
to the vortices JIB will induce a Lorentz force per unit 
length fi = (cj) /c)z x J. Thus we add to our model free 
energy of Eq. (fil) the term 



STr. 



J o j= i 



(21) 



Let us now estimate the effective glass exponent from 
the saddle-point configuration of Eq. (|l|), supplemented 
with Eq. (Eh) . A flux line will start to leave its columnar 
defect due to the current by detaching a segment of length 
z into the defect-free region, thereby forming a half-loop 
of characteristic transverse size r, see Fig. | (left). This 
costs the bending free energy efr 2 / ' z, and the lost pinning 
energy Vqz. This estimate uses on the average pinning 
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FIG. 4. Left: Two columnar defects of distance d with 
a vortex half-loop (solid line), and a double-kink excitation 
(dashed line), respectively. The Lorentz force Jl attempts 
to push the vortices to the right. Right: Double-superkink 
configuration, as required for variable-range hopping. The 
flux line "tongues" seek out a compatible low-energy, but not 
necessarily adjacent pin at distance R. 



in semiconductors fl The free energy cost of a double su- 
pcrkink of transverse size R and extension Z stems from 
the elastic term 2E k (d — > R) = 2E k R/d, and the dif- 
ference of pinning energies between the highest occupied 
state ei ss fj, and the empty site energy e 3 ; « /i + A(R) at 
distance R. In the presence of a small external current 
we can thus write the free energy aa3 



5T K w 2E k R/d + ZA{R) - f L RZ 



(24) 



The density of available states as a function of R with 
dimension D transverse to B (D = 2) on the one hand 



equals dP J 



H-A(B) 



g(e) de from considering the energet- 



ics, and on the other hand is given by « (d/R) D from 
simple statistics in real space. A(R) can then be deter- 
mined from the equation 



M+A(_R) 



g(e)de^R 



-D 



(25) 



potential Uq = (f/fe)- Taking the work by the Lorentz 
force into account, the free energy of the loop becomes 
approximatelyB 



&F>(r, z) w e x r 2 /z + U Q z - f L rz . 



(22) 



Optimizing this expression in the absence of external 
currents leads to the saddle-point configuration z* s» 
r * V^i/Uq, while for finite currents one finds r* s» 
cUo/(J(/>o), which after inserting into Eq. ( p2|) yields 



5J=Z 



3/2 



(23) 



We now identify the saddle-point free energy with 
the typical energy barrier U(J). For a sufficiently low 
current J\, the vortex half loop will extend to the 
nearest-neighbor pin of distance r* = d, and hence 
J i = cUo/(4>od). The flux line will then form a dou- 
ble kink of width w k = d^J e~i/Uo, which costs the free 
energy 2Ek = 2d-\/eiUo. Thus, for J\ < J < J c , with 
these expressions and Eq. (123) , we arrive at a current- 



voltage characteristics of the form (|19|) with the effective 
exponent p c g — 1. Recent experiments on BSCCO sam- 
ples have indeed found p c ff — 1 for intermediatep«irrents, 
supporting this half-loop nucleation picture.li3il3 



2. Variable-range hopping of flux lines 

For J < Ji, the most important thermally activated 
excitations will be double "superkinks" : The vortex will 
send out a tongue- like pair of kinks to a possibly very dis- 
tant neighbor, but with almost the same pinning energy; 
see Fig. || (right). Once the tongue is established, the su- 
perkinks move outward and the vortex hops to a new site. 
This process represents the flux-line analog of variable- 
range charge hopping transport between localized states 



We now proceed by optimizing Eq. (]2J) for J = 
first, which yields the optimal longitudinal extent Z* s» 
2Ek/[d (dA/dR)n*], while for a non-zero current one ar- 
rives at 



J^o/c w A(R*)/R* 



(26) 



which through inversion gives the typical hopping range 
R*(J). Inserting this last result back into Eq. (|24j),_we 
find the typical barrier for a jump, and consequentlyE3 

E = poJexpl-SJ^/T] = p Q Jcxp[-{2E k /T)(R*(J)/d)] . 

(27) 

This procedure can be followed numerically, starting 
with a given distribution of pinning energies g(e), as ob- 
tained from the energy-minimization algorithm. In order 
to determine the I-V exponent explicitly for a density of 
states of the form g(e) = j\e — /i| Soff , we insert this into 
Eq. (H) to get 



A(R) = 



and with Eq. 



Sett + 1 

7 



l/(s e ff + l) 



R -D/(s o{{ + 1) 



R*{J) = 



finally arrive at 

l/(D+s rff + l) 



7 



(28) 



*TJ ' (29) 



with p c ft given by Eq. (|2C|). Thus, the I-V characteristics 
will mainly be determined by the interaction-dependent 
effective exponent p c s, which itself depends on the effec- 
tive gap exponent s e g, i.e., the amount of depletion of 
the density of states in the vicinity of fi. With Eq. (|27j ) 
we can write the I-V characteristics finally in the form 



E = P oJcxp[-(2E k /T)(Jo/J) p ° !! } 



(30) 



with Jo = (c/0 o )7 1/(Sott+1) d 1/Poff . For the non- 
interacting case (with a constant density of states, i.e., 
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FIG. 5. Double-logarithmic plot (base 10) of the expo- 
nential factor R*(j)/d on SF>(j) vs the reduced current 
j = J(j>od/c, derived numerically from the site-energy dis- 
tribution g(e). Here, X/d — 1, and / = 0.8. 

Scff — 0), which is experimentally realized for very low 
magnetic fields / » 0.01, and thus X/ao — * 0, one obtains 
the Mott variable-range hopping exponent p c g = 1/3. 
Thus, we expect p e g to vary in the interval [1/3, 1] (the 
upper limit being given by s c ff - * oo), depending on the 
interaction range A, and the filling fractipjp—This has in- 
deed been found in recent experiments, EjO and partly 
in simulations. t2l 

With the distribution of site energies as obtained from 
our minimizing procedure, we can numerically evaluate 
the I-V characteristics for any form of g(e). This is done 
by integrating over g(e) starting at fi and using the re- 
lations (|5|) - (||) to obtain R*(J). A typical result 
derived from the distribution of pinning energies is de- 
picted in Fig. H which shows a double-logarithmic plot 
of the effective hop size R*(j)/d ~ 8F>(j) ~ U(j) as a 
function of the reduced current j — Jfad/c. According 
to Eq. (p9|), the slope, determined by linear regression, 
immediately yields the glass exponent p c s . We can thus 
extract p c s as a function of the filling fraction / and X/d 
from our data, and compare with our estimate of s e ff and 
the relation @. 

Apart from comparing the effective glass exponent 
measured in our simulations with experimental data, 
we may also predict certain features of low-temperature 
magnetization relaxation experiments that haxa.-been 
used to locate a possible Mott insulating phase.E3o As 
mentioned in Sec. |, these experiments identify a dip in 
the magnetization relaxation as the signature of a pos- 
sible Mott insulator. Assuming that the relevant relax- 
ation times are determined by variable-range hopping en- 
ergy barriers ln(i/t ) ~ U(J)/T = (2E k /T)(,Ul^ , 
one finds for the magnetization relaxation rate So'EjEj 

o_ dM _ dlnJ _ T 

dint dint p e sU(J) ' 1 ' 

and hence S ~ l/p e ft(/)- Previous studies of the magne- 
tization relaxation in the Mott insulator regime neglected 



this dependence of S on p c g, and rather investigated, 
quantum fluctuations on the Mott insulator transition,!!^ 
which may well play an important role in addition. Here, 
with our purely classical description, we argue that there 
may be some interesting physics in the behavior of p~g , 
which may already explain recent measurements of the 
relaxation rate near the expected Mott insulator tran- 
sition. Thus, from plotting p~~} vs /, we can at least 
qualitatively compare our results with these experimen- 
tal data. 



IV. DATA ANALYSIS 

In this section we present a detailed discussion of our 
simulation data. In the first part, we shall only use 
uniform, randomly distributed disorder, and discuss the 
properties of the Bose glass phase as a function of X/d 
and /. More specifically, we shall use X/d = 1/4, 1/2, 1, 5 
and 0.3 < / < 3. The main focus here lies on the 
question whether there exists a Mott insulating phase 
at / = 1, and if and when it is destroyed, as the inter- 
action strength and range between vortices is increased. 
We will also make contact between a number of experi- 
mental results and our simulations. In the second part, 
we shall in addition investigate a more correlated spatial 
disorder distribution, which only slightly differs from a 
triangular lattice. Here we want to investigate whether 
the Mott insulator can even exist at large ratios X/d. 

Before we start our analysis we would like to discuss 
finite-size effects. First, we want to point out that we 
do not perform a genuine finite-size scaling analysis of 
our data in this work. Rather, we intend to see quali- 
tatively and semi-quantitatively whether the Mott insu- 
lating phase may exist within the Bose glass phase, and 
if it is easily destabilized by strong long-range repulsive 
interactions. 

One may worry, however, about two issues in our sim- 
ulations: (i) There may remain artifacts from the under- 
lying grid structure, which we employ to make interstitial 
positions for the vortices available. We have tried to ren- 
der the grid finer by going from N = 1600 to N = 3600 
lattice points, keeping Nd = 100 fixed. However, we did 
not observe any significant lattice dependence in our mea- 
surements of \i or the total energy E in this case. Hence, 
we believe that our results should be largely insensitive 
to the underlying lattice representation and provide a fair 
approximation to a more realistic continuum description. 
We therefore used N — 1600 in most of our simulations 
presented here, (ii) One may worry that the number of 
defects is too small. To check for strong finite-size effects, 
we therefore chose N — 3600 and increased Nd from 36 
to 100 and 225. We did not detect any significant finite- 
size effects, as was also observed, in the simulations of a 
continuous defect distribution.E3 Also, we were able to 
make contact with Ref. [ [l2| by choosing / = 0.2 and 
/ = 0.4. We could reproduce the previous results in 
systems of similar sizes, e.g., the values of the chemical 
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FIG. 6. The filling fraction /*, at which 10 % of the vor- 
tices have left the defect positions, as a function of X/d. The 
resulting dashed line heuristically serves as a crossover line 
discriminating between the strongly pinned Bose glass (SBG), 
and the weakly pinned Bose glass (WBG) , respectively. 



potential fi and the effective exponent p c e . Hence, we be- 
lieve that for the purpose of our largely semi-quantitative 
analysis we can neglect artifacts stemming from too small 
lattice sizes or the underlying grid structure. 



A. Randomly Distributed Disorder 

1. Strongly vs weakly pinned Bose glass 

We start our data discussion by asking the question of 
how many flux lines become depinned as a result of the 
vortex interactions, depending on the filling fraction / 
and the interaction range X/d. As noted in the introduc- 
tion, within the Bose glass phase one may discriminate 
between regimes of strong and weak pinning (strongly 
and weakly pinned Bose glass, SBG and WBG, respec- 
tively), separated by a crossover line. The SBG is char- 
acterized by the fact that most vortices are pinned to 
columnar defects, i.e., single- vortex pinning is the domi- 
nant mechanism. In the WBG, on the other hand, only 
a fraction of lines is pinned by the defects, while many 
other vortices are held in place owing to the repulsive 
interaction exerted from other vortices. In this situation, 
one also speaks of vortex bundle pinning, i.e., an arrange- 
ment of flux lines bound to defects effectively form a cage 
for other vortices, and pin these through their mutual 
interactionBEsEZl For X/d <C 1 the crossovec-between the 
two regimes is expected to occur at B ~ B&B Once inter- 
actions become strong, however, the WBG should appear 
well below B$. 

Let us define the filling fraction /*, at which 10 % of 
the vortices are depinned in equilibrium as a result of 
the repulsive interactions. With /* we tentatively iden- 
tify the crossover line marking the change from SBG to 
WBG in our simulations, as a function of X/d. This line 
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FIG. 7. Distribution of (interacting) pinning energies g(e) 
vs e — n at / = 1 for four different interaction strengths. From 
top left to bottom right we have X/d = 0, 1/4, 1/2, 5. All plots 
include the data of 50 different realizations of the disorder. In 
the case of X/d — » 0, also the hard gap A separating occupied 
from unoccupied states, and the width of the random site 
energy distribution 2w are indicated. 



may thus be compared with the crossover line B*(T) as 
estimated in Rcf. and experimental results, see Refs. [ 
p^ , p^| . Of course, the specific criteria that are used to 
characterize the boundary between the SBG and WBG 
regimes are somewhat arbitrary, but the quantitative fea- 
tures of B*(T) should not depend too much on the exact 
manner this crossover line is defined. Figure ^ shows a 
plot of /* vs X/d. One observes a rather strong depen- 
dence of this line on X/d, which is shifted well below 
as soon as the interaction range reaches A ~ d. Only for 
X/d < 1/4 does the curve f*(X/d) remain above / = 1. 
This fact is already a first indication that the Mott in- 
sulator does not exist for X/d > 1/4, since it is defined 
by the localization of all vortices on the defects. This is, 
however, impossible, if already 10 % of the vortices are 
depinned as a result of their mutual repulsion at /* < 1. 
Notice that for X/d = 0, /* = 1.1 by definition. 



2. Site-energy distribution at the matching field 

In order to detect the appearance or disappearance 
of the Mott insulator more clearly, we now turn to the 
analysis of the (interacting) single-particle density of 
states g(e) at / = 1. Remember that by density of 
states we refer to the distribution of pinning energies in 
the flux-line picture, i.e., a histogram of the values of 
e, = n i n jV{ r ij) + U f° r k realizations of the disor- 

der distribution. For long-range interactions one expects 
a soft "Coulomb" gap in g(e) at the chemical potential fi, 
whereas for short-range interactions a hard gap A in the 
distribution of pinning energies should appear, indicating 
the presence of the Mott insulating phase. 

In Fig. |t] we show four different plots of g(e) vs e — \i for 
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different interaction regimes. Let us start at the top left 
of Fig. 0. In the limit X/d — > one expects a hard gap in 
the distribution of pinning energies, characteristic of the 
Mott insulator phase. This can be seen nicely in this plot, 
with a gap in <?(e) of width A = (Uk) — w separating the 
occupied states at e = — (Uk) ± w from the empty states 
at e = 0, with the chemical potential [i = —((Uk) — w)/2 
(by definition) located in the center. 

Let us now study whether the Mott insulator survives 
the introduction of fairly weak repulsive interactions. 
The smallest interaction range we can use on a lattice 
with N = leOCUand N D = 100, hence d = \J N/No = 4, 
is X/d = 1/4.E3 The top right of Fig. ^ depicts g(e) for 
X/d = 1/4. We observe a fairly flat density of states for 
the occupied sites, being somewhat depleted at e < [i 
and rises sharply for e > /i. However, the fact that all 
states below fi are continuously filled, implies that even 
rather short-range interactions suffice to destabilize the 
Mott insulator as a distinct thermodynamic phase, for its 
unique signature, namely the hard gap in the site-energy 
distribution, has vanished. Simultaneously, we observe 
that already 4 % of the vortices (in an average over 50 
realizations of the disorder) have left the defect sites, and 
moved to energetically more favorable interstitial sites. 
We have also looked at X/d 1/8 by increasing N to 
3600 and using only Nd = 56, and hence d « 8. Yet 
even in this case we found some states in the vicinity of 
fi and approximately 2 % of vortices on interstitial sites. 
This result strongly indicates that for the true Mott in- 
sulator to appear in the presence of a random, uniform 
disorder distribution, a necessary condition is A <C d. 
This incompressible phase may, therefore, never be truly 
realized in e xperiments, where usually A > d. We spec- 
ulate, though, that for a more correlated disorder dis- 
tribution the Mott insulator may persist to considerably 
larger interaction r anges ; this suggestion will be further 



investigated in Sec. |VB. 

We now turn to the two bottom graphs of Fig. |?]. For 
X/d = 1/2 (left) there is clearly no gap, but at e — /i ~ 
a narrow, soft gap begins to appear in the distribution 
of pinning energies g(e). In the bottom right figure for 
fairly long-range interactions X/d = 5, a wide soft gap of 
the form g(e) ~ |e — /i| Soff has fully developed. Such a 
Coulomb gap, though narrower, is also found for X/d = 1. 

We note that in this situation of strong repulsive in- 
teractions and large filling, the effect of the pinning sites 
themselves has become rather weak. The vortices are 
largely held in place by the cages formed by the other 
flux lines, and the overall structure should be rather un- 
stable against a collective drift motion, when an exter- 
nal current is applied. Quite obviously, we are now well 
in the weakly pinned Bose glass regime. In Fig. |^, we 
display a typical energy contour plot of an equilibrium 
configuration at / = 1 and X/d = 5. The top figure 
depicts the energy contour, where many circular lines 
indicate deep energetic minima, corresponding-to the de- 
fect/vortex configuration of the bottom plot.E3 For this 
case of fairly long-range interaction, the energetically fa- 
vorable sites include both the defect and high-symmetry 
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FIG. 8. Energy contour plot (top) for / = 1 and X/d — 5 
for the distribution of pinning sites (x) and the equilib- 
rium distribution of vortices (+)• Vortices occupying defects 
are marked by a star. Notice that even vortices that are 
not bound to defects are located at sites characterized by 
deep energy minima, see, e.g., the particles at coordinates 
(ar = 4,» = 4), (22,4) or (13,16). 



interstitial positions. We observe that pinning via the 
repulsive forces exerted by neighbors, which are still at- 
tached to defects, is very effective. As can be seen in 
Fig. ^|, the potential minima are roughly equally wide 
and deep, irrespective of the flux line occupying a defect 
or an interstitial site. 

In summary, by analyzing the density of states at / = 1 
we find that even for short-range interactions no distinc- 
tive Mott insulator phase appears. It becomes rapidly 
destabilized by the repulsive forces between flux linos, 
and only for X/d — > do we observe a finite hard gap in 
the density of states near fi, i.e., a Mott insulator. In the 
following, we shall investigate if the measurement of other 
observables corroborates this picture which we obtained 
by the analysis of the density of states. Furthermore, 
we discuss the qualititive explanation of recent experi- 
ments that identified (albeit rather broad) features near 
the matching field. 
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FIG. 9. Top: Bulk modulus en vs / on a log- linear plot. 
One observes a strong cusp at about / = 0.85 for X/d — 1/4 
(triangles). For X/d = 1/2 (open circles), the cusp is much less 
pronounced, and shifted to a lower filling fraction / m 0.55. 
Bottom: Reversible magnetization M vs / on a log-linear 
plot. For X/d = 1/4 there is a marked dip in the magneti- 
zation at / ps 1, whereas this structure is almost completely 
absent for X/d — 1/2, with only a shallow minimum at / « 0.6 
remaining. 



3. Bulk modulus and reversible magnetization 

Let us next turn to the determination of the bulk mod- 
ulus Cn as a function of the filling fraction /, while keep- 
ing X/d fixed. A divergence of en at / = 1 would indicate 
the appearance of the Mott insulator, as an incompress- 
ible phase (k = 1/cn = 0). In order to obtain en we have 
measured the chemical potential /x defined in Eq. (|l^) as 
a function of /, i.e., the H(B) curve. From this we obtain 
the bulk modulus by en = d/x/9/, the derivative being 
of course replaced by a discrete difference when analyzing 
our data. As we have seen in our study of the density of 
states, for A — * a hard gap of size A appears. In this 
case /x jumps from —A to as B is increased from 
to B§ plus one additional flux quantum. Thus, a diver- 
gence of en is trivially connected to a gap in the density 
of states. 

Fig. H (top) depicts the bulk modulus vs filling fraction 
on a log-linear plot. One observes that no divergence ap- 
pears at / = 1 for X/d = 1/4, yet a quite pronounced 
cusp occurs in en, shifted downwards to / ~ 0.85. For 
an even larger interaction range (X/d = 1/2) the cusp in 
/j, is displaced to even smaller / w 0.55 and less marked; 
for X > d, the bulk modulus becomes essentially constant 
over the entire range of /. In both cases depicted here, 
the location of the peak in en coincides with that value 
of / at which a sizeable number of vortices actually leaves 
the defect sites. The fact that the average bulk modulus 
is enhanced for larger A is simply an energetic effect, since 
a vortex entering the system has to overcome a higher en- 
ergy barrier due to the interactions. The behavior of en 
confirms again that no true Mott insulator exists even 



for comparatively short-range interactions X/d = 1/4 or 
1/2. We do observe, though, a distinctive "lock-in" struc- 
ture for X/d < 1, visible as a cusp in c\\. This feature, 
however, completely disappears for X/d> 1. 

In order to make contact with experiments, we mea- 
sured the total energy G, which yields the reversible 
magnetization via the thermodynamic relation M = 

-dG/dB N^dG/df. Fig. | (bottom) depicts M 

as a function of / on a log- linear plot for X/d = 1/4 
and X/d = 1/2, respectively. The data for X/d = 1/4 
show a pronounced minimum in M at / w 1, embed- 
ded in a slow logarithmic growth as / increases. The 
second plot for X/d = 1/2 hardly displays any structure 
aside from a shallow dip near / w 0.6. This feature is 
completely absent for X/d > 1, where only the In/ in- 
crease remains, resembling the magnetization curve of an 
unirradiated superconductor .EiJ'Ej'E3'E3 The observed be- 
havior of M at X/d =1/4 qualitatively agrees very well 
with recent measurements performed-by van der Beek et 
at in an irradiated BSCCO crystal,E3 where (at T m T\) 
a pronounced dip in M centered near B = B$ was found. 
The disappearance of this minimum in the experiments 
as T is increased towards the transition temperature may 
be at least partially due to the divergence of the London 
penetration depth A(T), which we capture in our data for 
large A, as well. In addition, the entropic renormaliza- 
tions studied in Ref. ^ will probably play a significant 
role. Our simulations cannot explain, however, the ma& 
netization minima found at / > 1 in BSCCO tapes.EJ 
We will return to this issue in Sec. [v| 

We would like to point out that the maximum in en 
occurs at lower value of / than the minimum in M (com- 
pare, e.g., our data for X/ d = 1 /4), which can be under- 
stood as follows (see Sec. Ill D| ) . Starting from the rela- 
tion M = B — H and n ~ H — H Cl , the bulk modulus 
may be written as en ~ dH/dB = 1 — dAI/dB; thus 
the maximum of en occurs at the location of the steep- 
est negative slope in M(B), which has to be at a smaller 
B than the minimum in M. Only in the "true" Mott 
insulator phase will the singularities in both en and M 
coincide precisely at B = B$. 

In summary, the measurements of the bulk modulus 
and the reversible magnetization confirm the scenario 
that no Mott insulator phase exists at finite X/d. There 
are, however, interesting "lock-in" structures present at 
X/d < 1, which qualitatively agree with recent experi- 
mental findings. We interpret these as remnants of the 
Mott insulator phase, but not as true signatures of this 
distinctive thermodynamic phase itself. 



4- I-V characteristics and magnetization relaxation 

In the last section of our data analysis we now turn 
to the transport properties of our model, which can be 
evaluated from t he sin gle-particle density of states as de- 
scribed in Sec. 



HIE. The calculation of the effective 



transport exponent p e ff by such means was first used 
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in simulations in Ref. [ [T^], where small filling frac- 
tions only were analyzed and thus no "interstitial" grid 
was used as in the present study, but rather a continu- 
ous distribution of defect sites. Thus, each vortex was 
always pinned by an energetically favorable defect site 
(single-vortex pinning) and hence the equilibrium con- 
figurations should be stable against thermal fluctuations 
or small currents (strongly pinned Bose glass regime). 
In our case, where we want to analyze / > 1 as well, 
we need to see whether our equilibrium configurations 
are sufficiently stable against such perturbations. But as 
the above energy-contour plot (Fig. |J) demonstrates, the 
potential minima for the vortices both at defect and at 
interstitial positions are about equally deep, and gener- 
ally well-separated from each other. This gives us con- 
fidence that the single-particle density of states may in- 
deed be used to infer low-current transport properties 
in the variable-range hopping regime. It also justifies a 
posteriori the quality of our grid description of the prob- 
lem: The most important low-energy interstitial sites are 
correctly taken into account. 

We have determined the effective glass exponent de- 
fined by the equation 



E = p Jexp[-(U c /T)(J 1 /J)*«] , 



(32) 



the 



via employing the method described in Sec. [HE 
range 0.5 < / < 3 for X/d = 1/2, 1, 5. We were not able 
to extract p e g from the density of states for X/d = 1/4, 
since the absence of even a small Coulomb gap (cf. Fig. |t], 
top right) renders the definition of s c g impossible. Our 
findings are that for interaction ranges X/d — 1,5 we ob- 
serve a slow, monotonous rise in p e g as / is increased: For 
X/d = 1 from p eS (f = 0.5) « 0.54 to p cS (f = 3) » 0.74, 
and for X/d = 5 fromp e ff(/ = 0.5) « 0.7 to Pes(f = 3) w 
0.8. It is interesting to compare these results with the 
earlier simulations at low filling fractions£3 For instance, 
there at / = 0.1 and X/d — 1 an effective exponent of 
p cff w 0.5 was found, which then decreased almost to the 
non-interacting Mott value p c g = 1/3, as / was increased 
towards 0.8. This, however, was due to the fact that the 
system had to accommodate with the underlying ran- 
domness, since no favorable interstitial sites were avail- 
able. In our case, with these interstitial sites being acces- 
sible, the disorder effects are screened by the interactions 
as the vortex density is increased, which leads to stronger 
correlations, and presumably more mean-field like behav- 
ior (p c ff — ► 1), because of the more uniform spatial distri- 
bution of flux lines. Thus, our simulations are presum- 
ably more realistic for higher filling fractions, and the 
tendency of increasing transport exponents with larger 
magnetic fields is experimentally confirmed in recent I-V 
measurements on YBCO.t3 In these experiments, a con- 
tinuous increase of p e g with B from p c s(f ~ 10~ 2 ) = 1/3 
to Pefi(f ~ 1) ~ 1 was found. 

We plot our data for p~^ in Fig. @. For X/d = 1/2 
we find a marked minimum in p~ s at / w 1.4, which is 
presumably due to a delicate interplay of disorder and 
interaction effects. Assuming that relaxation times are 




FIG. 10. Inverse transport exponent l/p e ff vs / for 
X/d — 1/2, 1,5. For \/d — 1/2, the maximum in p e ff is lo- 
cated at B m I.4.B0 . 



determined by the variable-range hopping energy barri- 
ers, we found in Sec. HIE that the magnetization relax- 
p~fi(f). Thus, the minimum we find in 



ation rate S 

p~g for X/d = 1/2 might explain low- temperature re- 
laxation experiments on YBCO samples,LJ as well as on 
Tl-compounds,E.a for which a marked minimum in S(B) 
at B w 1.4£>$ wa reported. With the knowledge we have 
gathered from the above measurements, we would argue 
that the minimum detected in the magnetization relax- 
ation experiments is not really a signature of a true Mott 
insulating phase, but rather a remnant of the Mott insu- 
lator in the weakly pinned Bose glass near / = 1. 

In summary, our simulation data suggest that in a 
sample with a uniform random distribution of colum- 
nar defects, a Mott insulating phase is not observable 
as a consequence of the repulsive interactions between 
vortices. It may, however, persist up to considerably 
larger interaction ranges, if the disorder distribution is 
chosen in a more correlated manner, e.g., almost match- 
ing the triangular Abrikosov lattice of a pure type-II 
superconductor £3 This way the strong repulsion of vor- 
tices on defects, which are located very closely to each 
other, is avoided, and the interactions basically lead to 
an overall constant upwards shift of all the site energies. 
We are going to investigate this possibility in the follow- 
ing section. 



B. Spatially Correlated Disorder Distribution 

In order to investigate the presence of a Mott insu- 
lator phase up to larger interaction scales in a spatially 
more correlated disorder potential, we have constructed a 
trial disorder distribution as follows: The triangular lat- 
tice points with mutual distance 4a are looked up on our 
underlying grid, and then, with equal probability, the 
defect is either put on the six nearest neighbors of the 
original point, or it remains on the original site. Thus, 
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disorder only comes in as a slight distortion of a per- 
fect triangular defect lattice. Related simulations using 
ideal triangular and square disorder lattices were per- 
formed by Anghelache et al.jEi with the remarkable re- 
sult that at some values of / and X/d, a quadratic defect 
array is apparently more effective in pinning the vortices 
than a triangular lattice of pins. We remark, however, 
that the present Monte Carlo algorithm was specifically 
designed for disordered systems, and does not always 
faithfully reproduce the correct ground state in ordered 
structures (typically, regularly ordered subdomains re- 
main separated by dislocations, which cannot be removed 
by the merely local processes allowed for by the algo- 
rithm). Such regular lattices of pins have recently been 
produced in high-T c superconductors, and also been in- 
vestigated experimentally, where the main focus lay an 
the enhancement of the irreversibility line, i.e., T c (B)x3 
As for temperatures at least, the pinning efficiency of 
such regular pin arrays is much improved, we expect the 
crossover line between the weakly and strongly pinned 
Bose glass regimes to be drastically shifted upwards. We 
also remark that the non-equilibrium properties of such 
regular arrays under conditions of strong drive were stud 
ied numerically by Reichhardt et alrH 



1. Strongly vs weakly pinned Bose glass 

We first ask again the question, how many vortices 
have left their defect positions as a function of X/d? It 
is interesting to note that by looking at the defect oc- 
cupancy at / = 1, we find that even for X/d = 5 only 
roughly 1 % of the vortices have left the defect posi- 
tions, preferring an interstitial position. At X/d < 1 the 
vortices are always located on defect sites. Thus, the 
crossover line between the strong and the weak Bose g lass 
which we have defined pragmatically in Sec. IV A 1, re- 



mains above / = 1 for the entire range of interactions we 
study here. This is in striking contrast to Fig. ^, where 
we observed a quickly decreasing crossover line as soon 
as X/d was increased above 1/4. This indicates that a 
more correlated type of disorder enhances single-vortex 
pinning drastically, which has already been verified in 
experiments.Ej Moreover, it renders the occurrence of a 
Mott insulating phase for X/d < 1 possible. 



2. Site-energy distribution at the matching field 

As in the previous section, we set out to detect a pos- 
sible Mott insulator by looking for a hard gap A near the 
chemical potential in the density of states 17(e) (distribu- 
tion of pinning energies). Figure [II] shows four different 
plots of the density of states at / = 1 for our more cor- 
related disorder distribution. Starting with the top left 
of Fig. m, we see that at X/d = 1/4 there clearly ex- 
ists a broad hard gap in 17(e) near the chemical potential 
\i. All vortices are localized on defects in all 50 realiza- 
tions of the disorder and thus a true Mott insulator phase 
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FIG. 11. Distribution of pinning energies g(e) vs e — /j, for 
four different interaction strengths, using the correlated dis- 
order distribution. From top left to bottom right we have 
X/d= 1/4,1/2,1,5. 



emerges in this interaction regime. This is in contrast to 
the uniform random disorder distribution of the last sec- 
tion, where for X/d = 1/4 and even smaller values the 
Mott insulating phase was already destroyed (cf. Fig. [7], 
top right). 

The top right figure with X/d = 1/2 shows a signifi- 
cantly smaller, but still finite gap between the occupied 
low-lying states and the vacant high-energy states. Only 
at X/d = 1 (left bottom) has the gap closed and a soft 
Coulomb gap structure has emerged instead. However, 
all vortices are still bound to defects. Only for X/d = 2 
(not depicted here) do the first vortices leave their defect 
positions in favor of low-cost interstitial sites. Finally, for 
X/d = 5 the Coulomb gap structure only insignificantly 
differs from the one found for the uniform random dis- 
order distribution (see Fig. (7] bottom right); the disorder 
is now effectively screened by the long-range interaction 
between vortices. 

We thus conclude that the Mott insulator phase can 
survive the introduction of inter- vortex interactions if the 
defect structure has a more correlated form, as studied 
in our example. We find that the Mott insulator can still 
be found up to X/d ~ 1, in sharp contrast to the com- 
pletely random defect distribution, in which case even for 
X/d ~ 1/8 no incompressible phase could be identified. 
Let us now see how the Mott insulator phase appears in 
the measurement of the bulk modulus and the reversible 
magnetization. 



3. Bulk modulus and reversible magnetization 

In order to see the effects of the Mott insulator on the 
bulk modulus and the reversible magnetization, we con- 
sider the interaction regime X/d — 1/2, where we still 
observe a finite gap in the density of states. In Fig. |lj 
we show the chemical potential /1 vs /. This function 
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FIG. 12. Chemical potential n (top left), bulk modulus en 
(top right), total energy G (bottom left), and magnetization 
M (bottom right) vs filling fraction / on log-linear scales. The 
interaction length scale for all plots is X/d = 1/2. 



sulator is visible. X/d — 1 roughly represents the border 
line between the occurence of a hard gap in the density 
of states and the soft Coulomb gap. The cusp in c\\ near 
/ = 1 is still rather narrow, but much less pronounced, 
and also the dip in the magnetization is less strong. For 
X/d = 5, the bulk modulus remains flat over the entire 
range of the filling fraction, and the magnetization shows 
no significant dip. 

In summary, we have seen that a more correlated dis- 
order distribution enhances single-vortex pinning drasti- 
cally, and may thus even show a Mott insulating transi- 
tion for moderate interaction ranges. The Mott insula- 
tor phase should be observable in measurements of the 
reversible magnetization, with a sharp dip precisely at 
/ = 1 being its clear signature. However, this phase 
might well be rather unstable with respect to thermal 
fluctuations, because correlated vortex hops would be fa- 
cilitated by the more regular defect distribution. 

V. CONCLUSION 



displays a distinct jump near / = 1, slightly smeared 
out due to finite-size effects. Hence, the bulk modu- 
lus en = dj-i/df shows a strong narrow maximum near 
/ = 1, which in the limit of an infinitely large system 
would presumably yield a ^-function peak. Notice that 
the bulk modulus stays fairly constant below and above 
/ = 1. This behavior in the bulk modulus needs to be 
contrasted with that of en for the uniform random dis- 
order distribution and X/d — 1/4 (cf. Fig. ^top). There, 
at an even shorter interaction scale, only a broad and 
much flatter cusp in en was found, and its maximum 
was even shifted to a lower filling fraction of / « 0.85. 
This had indicated the absence of a Mott insulator phase, 
although the system near / = 0.85 was clearly fairly stiff 
with respect to introducing new vortices. Here, however, 
the narrow and sharp cusp in en and the gap in the den- 
sity of states strongly indicate the presence of the Mott 
insulator. 

The bottom left plot of Fig. [l^ shows the total en- 
ergy G vs /. Clearly, there is a marked change in the 
slope of G as the filling fraction is tuned through / = 1. 
Consequently, there is a strong dip in the magnetization 
M = —dG/dB at / = 1, as can be seen in the bottom 
right plot. This minimum is much more pronounced and 
sharper as compared to those found for X/d = 1/4 or 1/2 
in the uniform random disorder case (compare the broad 
dips in Fig. || bottom). Also note that both the peak in 
en and the dip in M occur at the same location, namely 
precisely at / = 1, as is indicative of th e true Mott in- 
sulator phase (see our discussion in Sec. pIPp . To our 
knowledge, such a sharp dip indicating a real Mott in- 
sulating transition has not been seen in experiments so 
far. 

Let us just briefly describe the other interaction 
regimes. For X/d = 1/4 of course the "divergences" in 
c\\ and M at / = 1 feature even stronger than in the case 
depicted in Fig. O and, therefore, also here the Mott in- 



The aim of our study was to investigate how the re- 
pulsive forces between vortices affect the properties of 
the Bose glass phase, and especially the proposed exis- 
tence of a Mott insulator phase, as the magnetic field is 
varied near the matching field B<$. Specifically, we set 
out to explain recent experiments which apparently have 
found evidence for the Mott insulator transition in irra- 
diated high-T c -superconductors. Our main findings for a 
uniform random defect distribution, which is realized in 
most experiments, are: 

• At the matching field B$, the Mott insulator, be- 
ing characterized by a hard gap in the distribu- 
tion of pinning energies g(e) near the chemical po- 
tential /i, only exists for extremely short-range in- 
teractions X/d — ► 0. With increasing X/d, the 
gap quickly fills and, therefore, the Mott insulator- 
phase as well as the singularities of the accompany- 
ing phase transition are destroyed. For long-range 
interactions A ^ d, a soft Coulomb gap described 
by g(e) cx \e — /i| Soff emerges with an effective gap 
exponent s c g > 0. 

• By measuring the reversible magnetization and the 
effective I-V exponent p e g, which we argue is re- 
lated to the magnetization relaxation rate—we-are 
able to explain recent experimental resultsE3~L.il as 
remnants of the Mott insulator for A < d, but not 
as true signatures of this thermodynamic phase it- 
self. More specifically, the dip in the magnetization 
at / ^ 1 and the minimum in the relaxation rate 
S near B sa lAB^ found in the experiments can 
be reproduced in our simulations, albeit with a si- 
multaneous absence of a hard gap in the density of 
states. 

• By tuning / in the interval [0.5,3] for X/d > 1 we 
find a slowly increasing effective glass exponent p e ff 
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with /, which is in accord with recent measurement 
on YBCO films.Pil 

• We can identify a crossover line discriminating be- 
tween the strongly and weakly pinned Bose glass 
regimes, respectively. This curve is close to the 
matching field for X/d < 1/4, but decreases rapidly 
as the interaction range is increased. 

Our simulations cannot explain the experimental find- 
ing of Li et aZ.,Eil who observed a minimum in the re- 
versible magnetization at B w 1.6_B$ in two samples of 
irradiated BSCCO. It is, however, intrinsically difficult 
to guess the correct value of -B$ , which is estimated from 
knowing the average flux density of the heavy-ion beam 
and some guesses of how effectively these ion columns can 
create pinning tracks for the vortices. Slight deviations 
of the estimated from the true defect density may there- 
fore be responsible for this apparent-discrepancy of the 
minimal magnetization and their _B$.Eil In this case there 
would be no contradiction to ouwjesult and that of the 
experiments of van der Beek et al.c3 Of course, we cannot 
exclude that thermal fluctuations might also play a role, 
as these were completely neglected in our algorithm, but 
these would rather be expected to further smoothen any 
structures in the magnetization. 

Can it be understood microscopically why the Mott 
insulator is destroyed for even short-range interactions? 
Due to the random distribution of pinning sites, there 
is a small probability in the Poisson distribution that 
two, three, or even four defect sites are clustered closely 
together, say, three defects on a triangle of our grid. 
Although these defects present favorable pinning sites 
for vortices, their mutual (though short-range) repulsion 
may lead to the effect that one vortex actually leaves 
in favor of a high-symmetry interstitial site, while the 
others remain on the pinning centers. Thus, even for 
short-range interactions, a small but finite percentage of 
flux may never sit on defect sites due to the clustering of 
defects. Consequently, the Mott insulator, in which all 
defects are occupied by vortices, is destroyed. 

To cancel this effect we also used a trial disorder dis- 
tribution, which is basically a slightly distorted triangu- 
lar defect lattice. We presumed that in this case single- 
vortex pinning should be enhanced and thus the Mott 
insulator should be visible for larger interaction regimes. 
Here our results are as follows: 

• The crossover line discriminating between the 
weakly and the strongly pinned Bose glass is drasti- 
cally shifted upwards. It remains above / = 1 even 
for X/d = 5. 

• At the matching field the Mott insulator can be 
found up to X/d < 1/2, marked by hard gap in 
the density of states near the chemical potential. 
The gap fills, and a soft Coulomb gap emerges for 
X/d> 1. 

• For A < d we find a sharp and narrow peak in 
the bulk modulus at / = 1; simultaneously, we ob- 



serve a distinct sharp minimum in the magnetiza- 
tion at B = _B$. Both signatures are much more 
pronounced than even for X/d = 1/4 in the uniform 
random disorder case. The fact that both singulari- 
ties occur simultaneously at / = 1 is another strong 
indication of the existence of a true Mott insulator. 

Our simulation results indicate that in most experi- 
mental samples studied so far, which have hardly any 
correlations in their disorder distribution, the Mott in- 
sulator is presumably never realized. In most high-T c - 
superconductors the interaction length scale A far below 
T c is of the order 1000 A. The experiments of Li et al. and 
van der Beek et al. used a fairly high irradiation density, 
such that the average defect distance was of the order 500 
A, and thus X/d w 2. Presumably a Mott insulator could 
be produced with a far weaker radiation dose such that 
d^$> X. Another possibility would be using a more corre- 
lated disorder distribution, like in our simulations, such 
that repulsive forces from very close vortices occupying 
defects can be avoided. It would be interesting to investi- 
gate both possibilities experimentally. Furthermore, the 
inclusion of thermal fluctuations in the investigation of 
the Bose glass properties would be highly desirable. How- 
ever, this requires highly efficient algorithms that are able 
to deal with three-dimensional vortex line fluctuations in 
a disordered environment. 
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